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Abstract 

We propose a growing network model that consists of two tunable mechanisms: 
growth by merging modules which are represented as complete graphs and a fitness- 
driven preferential attachment. Our model exhibits the three prominent statistical 
properties are widely shared in real biological networks, for example gene regulatory, 
protein-protein interaction, and metabolic networks. They retain three power law 
relationships, such as the power laws of degree distribution, clustering spectrum, 
and degree-degree correlation corresponding to scale-free connectivity, hierarchical 
modularity, and disassortativity, respectively. After making comparisons of these 
properties between model networks and biological networks, we confirmed that our 
model has inference potential for evolutionary processes of biological networks. 

Key words: power law; biological networks; scale free; hierarchical modularity; 
disassortativity 



1 Introduction 



Structural data of large-scale biological networks such as gene regulatory net- 
works [T1I21I3] . protein-protein interaction networks FilloIRi 7 . and metabolic net- 
works [8f9] have been accumulated by the active investigations in recent years, 
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and heterogeneous connectivity, which means that a few nodes integrate a 
great number of nodes and most of the remaining nodes do not, are ubiqui- 
tously found as a striking property of the structures Pp|7f5ll^lll0flllll2lll3fl4j . 
Heterogeneity exhibits the scale-free feature [TofTo] . and differs from homoge- 
nous connectivity, which means that most nodes have the same degree (number 
of edges), is observed in classical network models such as regular lattices and 
the Erdos-Renyi random network [T7] . The scale-free feature is in degree distri- 
bution, defined as the probability of the existence nodes with degree k, which 
follows the power law, P(k) oc k' 1 with 2 < 7 < 3. Since the scale-free fea- 
ture is independent of living species [8"fll|ll2|ll8j , growth mechanism which can 
generate such statistical properties are paid much attention, and are expected 
to elucidate the evolutionary and growing processes of the networks. 



In previous works, the Barabasi- Albert (BA) model [16J is well-known as a 
model for scale-free networks, and is frequently embedded among variety of 
mathematical models [19] . Using the growth mechanism by adding new nodes 
with preferential attachment (PA) which means that high-degree nodes get an 
even better chance to attract new edges, the BA model generates scale-free 
networks which indicate power-law degree distributions with the fixed degree 
exponent —3; thus P(k) oc k~ 3 [T(3lll8|l2U] . However, the recent detailed statis- 
tical analyses of biological networks manifest some discrepancies between the 
networks generated by the BA model and biological networks. In particular, 
there are two striking properties as follows. 

One of the properties is a small- world feature [2T]. This feature is reflected 
in high clustering coefficients C, which denote the density of edges between 
neighbors of a given node, and implies modularity of networks [22)123] . The 
modular structures are actively investigated with statistical approaches, and 
it is found that clustering spectra, defined as correlations between degree k 
of a given node and the clustering coefficient C of the node, follow the power 
law with exponent —1; thus C(k) oc k~ x [T3"|14f24j . The power-law spectrum, 
observed in biological networks, suggests a hierarchical structure of modules 
[23]. However, BA networks do not exhibit a power-law clustering spectrum, 
indicating the absence of the hierarchical structures [25T2 6J. In order to fill 
the gap, we propose a model (hereinafter called MM model) with growth 
mechanism by merging modules and PA from the BA model. The MM model 
shows power-law clustering spectra and power-law degree distributions with 
arbitrary degree exponents, demonstrating that our model can reproduce finer 
details of biological networks 



The other of the properties is a heterogeneous degree-degree correlation 
defined as correlation between degree A; of a given node and average degree k nn 
of neighbors of the node. In biological networks, the degree-degree correlations 
follow power law; k nn (k) oc k v with a negative exponent, — 1 < v < 0, em- 
pirically found [T31ll4f29] . In general, the properties are called disassortatvity 
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|30j . It is reported that the BA and MM model have no disassortativity [2~Tf2"8] . 
Pastor-Satorras et al. [28] and Barrat et al. [31] respectively argue that dis- 
assortativity emerges with competitive dynamics [32] and weight- driven dy- 
namics [33], suggesting that growing mechanism needs to include fitness such 
as spatial distance [34J, aging [35J, and so on. 

In this paper, we propose a model which maintains the following statistical 
properties: the scale-free feature, the hierarchical modularity, and the disas- 
sortativity. For the proposal, we start with a modification of the MM model 
that has the scale- free feature and the hierarchical modularity [27]. We embed 
a tunable fitness-driven (FD) mechanism into the MM model, and provide 
the numerical and analytical solutions for the statistical properties of model 
networks. 



2 Model 



Here we present an evolving network model. Our model includes growth by 
merging modules and the FD mechanism [see Section [TU] for biological impli- 
cations of the model]. 

Growth by merging modules — A network grows through joining new modules 
to existing nodes of a network over time [see Fig. [1] (C)]. m/a is the first 
control parameter of the model where a and m denote the number of nodes 
of the complete graph and the number of merged node(s), respectively. Please 
note that the process develops without adding extra edges 



Fitness- driven (FD) preferential attachment (PA) — The standard PA mech- 
anism of the BA model is the probability that node % is chosen to get an edge, 
and is proportional to the degree of node i; hence Ilf A = ki/J2jkj, where 
ki is the degree of node i. The mechanism only considers the degrees at the 
nodes. Here, we additionally consider the probability that node % is selected 
according to degree k^ and fitness / i; and express the probability as 

n = k + fi (1) 



Updating rule of fitness — Moreover, we consider the change of the fitness of 
node i. When node % is selected by using the FD-PA mechanism given in Eq. 
(CQ), the fitness of node % increases as follows: 

fi^fi + Z, (2) 
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(A) (B) (C) 

Fig. 1. Schematic diagram of growth process of our model network with a = 3, 
m = 1, and £ = 5. (A) Selection of node(s) by the FD-PA mechanism, Eq. (fT]). The 
filled node is selected by the attachment. Each number in the figure indicates the 
fitness of a corresponding node. (B) Updating of the fitness. The selected node's 
fitness increases according to the updating rule, Eq. ([2]). (C) Merging new modules. 
As a result, (a — m) new node(s) are added, filled with gray, and take zero as their 
initial fitness. 

where £ takes a constant value, and is the second control parameter in the 
model, indicating the strength of incidence of the fitness. 

Taken together, our model networks are generated by the following procedures. 

i) We start from a module that is a complete graph consisting of a (> 3) 
nodes. For fitness, we assign zero to all nodes in the module. 

ii) At every time step a new module with the same size is joined by merging 
to an existing m (< a) nodes. 

iii) When merging the module, the FD-PA mechanism, Eq ([1]), is used to 
select m old nodes [see Fig. [I] (A)]. Then, the fitnesses of the old nodes 
are increased using the updating rule, Eq. ([2]) [see Fig. [I] (B)]. At the 
end, zero is assigned to the fitnesses of the new nodes [see Fig. CD (C)]. 
Resultant duplicated edges between the merged nodes are counted and 
contribute to the FD-PA mechanism in the succeeding steps. 

When a > m and £ > 0, the model network evolves in time steps. In addition, 
our model is equivalent to the BA model with the specific condition: a = 2, 
m — 1, and £ = 0. 



3 Degree distribution 

In this section, we present the analytical and the numerical solutions for the 
degree distribution of our model. The degree distribution is an important 
statistical property for the characterization of the networks, and denotes the 
existence probability of nodes with degree k. The degree distribution [T8f2"0] 
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is denned as 



1 N 

^(*) (3) 

i=l 



where S(x) and N are Kronecker's delta function and the total number of the 
nodes, respectively. 



3.1 Analytical solution 



In order to describe the degree distribution, we take the continuous mean-field 
approach [16.26J. The standard approach can not be applied directly due to 
the inclusion of the fitness updating. We express the degree and the fitness as 

Fi = ki + f u (4) 



and investigate a time evolution of Fi. Since the fitness updating rule is rep- 
resented as fi <— fi + £, Fi satisfies 

F <=(j^T + 1 )*'-f' (5) 



indicating the proportional relationship between Fi and fcj. 
The time evolution of Fi is described as 

f = m(a -l + f)^, (6) 

where ~ [a(a— 1) + m£]i. The solution of the equation with Fi(t = s) = 

(2) + £ = ^(°>0 as an initial condition for Eq. ([6]) is 

F^)=A(a,0(^, (7) 

where (3 = [m(a — 1 +^)]/[o(o — 1) + m£]. Because s/t denotes the probability 
that i 7 ^ is larger than a given F, Equation (J7J) is rewritten as 

P(>F) = A(a,£) 1/p F- 1 / f) , (8) 
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which corresponds to the cumulative probability. From Eq. (jSj), the probability 
distribution for F is given as 

P(F) = - A P(> F) = ^jt!l F -^ (9) 



where 7 = + 1. Finally, substituting Eq. (jSJ) into Eq. (Q, we get the 

degree distrubution 

P(k)~B{a,Z,>y)k-\ (10) 



where B(a, £, 7) = A(a, £) 7 1 [^/(a — 1) + 1] 7 / (7 — 1). The degree distribution 
obeys the power-law with the degree exponent 

(q + m)(q- 1) + 2m£ 
7 ~ m(a-l+0 ' 1 j 



demonstrating that our model network represents the scale-free feature. 



3.2 Numerical solution 



In order to confirm the analytical predictions, we performed numerical sim- 
ulations of networks generated by using our model described in Sec. [2j We 
show degree distributions with a = 5, N = 50000, different m G [1,4] and 
£ G [0, 15] in Fig. El The degree distributions follow the power law, reflecting 
the scale-free feature. 

Figure [2] (A) shows the degree distributions with m = 2 and the different £. 
The degree exponents decay with increasing £. Figure [2] (B) shows the degree 
distributions with £ = 7 and the different m. The degree exponents decay with 
increasing m. We show excellent agreement between the numerical results and 
the theoretical predictions, demonstrating that 7 is function of m/a and £. 



4 Degree entropy 



As shown in Eq. (fTTj) and Fig. [2j the degree exponent 7 is variant with £ 
and m/a. We consider the degree entropy [36] to characterize parameter de- 
pendency for the heterogeneous connectivity of our model. Assuming that 
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Fig. 2. Degree distributions P (k) with a = 5 and N = 50000 (shifted for clarity). 
Different symbols correspond to the different numerical results. Solid lines show 
relationship oc /c~ 7 , where 7 is predicted by Eq. (fTTj) . (A) £ dependency with m = 2. 
(B) m/o dependency with £ = 7. 

X^a^Ti 1 -P(^) = 1, the degree entropy is denned as 

N-l 

H = - J2 P(k)\nP(k), (12) 
fc=i 



and provides an average measure of a network's heterogeneity because it relies 
on the diversity of the degree distribution. The increment of degree entropy H 
corresponds to the increment of the heterogeneity for complex networks [36] . 

Figure [3] shows the numerical results of the degree entropy of the model net- 
work with a = 5, N = 12800, different m G [1,4], and different £ e [0, 15]. 

In Fig [3] (A), we show the degree entropy H with the different m/a and 
£ in three dimensions. The degree entropy H increases for the larger m/a 
and the smaller £, and decreases in the opposite direction. In Fig [3] (B), the 
projections of the plot to the (m/a)-H plane of Fig. [3] (A) are shown. The 
correlation between H and m/a becomes strong with increasing £. In Fig [3] 
(C), the projections of the plot to the £-H plane of Fig. [3] (A) are shown. The 
degree entropy decays with increasing £ in the case of constant m/a. 



5 Clustering spectrum 



The clustering spectrum is a well-known statistical property which reflects the 
hierarchical modularity of networks [T8"1I2T)] . and is defined as 
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Fig. 3. Degree entropy H with a = 5 and ./V = 12800 (Axis direction changed for 
clarity). (A) 3D plot of the entropy versus the variables £ and m/a. Projections of 
the plot to the (m/a)-H plane (B) and the £-H plane (C). 

where S(x) denotes Kronecker's delta function, and Cj corresponds to the 
clustering coefficient [2T], defined as 



and represents the density of edges among neighbors of a given node, where 
hi and M\ denote the degree of node % and the number of edges among the 
neighbors, respectively. 

First, we give the analytical solution for the clustering spectrum of our model. 
Since our model grows through merging modules (see Sec. [2] and Fig. [[]), the 
number of edges of node i among neighboring nodes is approximately described 
as [27] 



Mi 



2Mi 



(14) 






(15) 



where Si corresponds to the number of the selection of node i with the PA, 
as in Eq (CQ). In our model, because the degree of node i is expressed as 




k k 



Fig. 4. Cluster spectra C(k) with a = 5 and N = 50000. Different symbols cor- 
respond to the numerical results. Solid lines show the relationship oc k^ 1 . (A) £ 
dependency with m = 2. (B) m dependency with £ = 7. 

fcj = 5i(a — 1), Eq. (fT5|) is rewritten as 

M, ~ ^fc, (16) 



indicating the proportional relationship between and fcj. Finally, substitut- 
ing Eq. ffTBT) into Eq. flHJ), we get the clustering spectrum 

C(jfe) ~ oc k-\ (17) 



The clustering spectrum follows the power law with the exponent —1 when 
a > 3, reflecting the hierarchical modularity of our model network. 

Next, we present the numerical results of our model for the clustering spectrum 
to verify the analytical solution. In Fig. HI we show clustering spectra with 
a = 5, N = 50000, different m G [1,4], and different f G [0, 15]. 

Figure 0] (A) shows the clustering spectra with m = 2 and the different £. 
The spectra follow power law with the exponent —1 despite changing £. Fig- 
ure H] (B) shows the clustering spectra with £ = 7 and the different m. The 
spectra follows power law with the exponent —1 for large k. For small degree 
k, the cut-off is prominent with increasing m due to running away from the 
approximation of Eq. (TT51) . 

The clustering spectra follow the power law with the exponent more or less 
equal to —1, reflecting the hierarchical modularity of networks. 
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6 Global clustering coefficient 

As seen in Fig. [3], degree entropy H is sensitive to £ and m/a. To examine 
parameter-dependency of the modularity (clustering property) for the whole 
network, we utilize a global clustering coefficient C. The coefficient is defined 
as [21] 

1 N 

C=^- T Y J C i . (18) 

A higher value of C means that the network has higher modularity. 

Figure [5] shows the numerical result of C in our model network with a = 5, 
N = 12800, different m G [1,4], and different £ G [0, 15]. 

In Fig. [5] (A), we show the C with different m/a and £ in three dimensions. C 
tends to increase more for smaller m/a and larger £. 

Since the small or large m/a give rise to different effects on C, Figure [5] (A) 
has a valley in the middle of m/a. In the case of small m/a, C can remain high 
because the modules combine via few common nodes. In the case of large m/a, 
in contrast, C leniently decreases because the networks tend to be randomized 
since the modules combine via many common nodes. In this case, however, FD 
mechanism helps to increase the fitness of the nodes in the modules, inducing 
a formation of a cluster with high-edge density. As a result, C increases with 
m/a. Due to the trade-off mechanism, the biphasic graph is shown in Fig. [5] 
(B). In addition, C is less sensitive for large £. In Fig. [5] (C), we show the 
projection of the plot to the £-C plane of Fig. [5] (A). 



7 Degree-degree correlation 

The degree-degree correlation is a statistical property that characterizes as- 
sortativity of the networks, and represents the average degree of neighbors of 
nodes with degree k [181120] . The correlation is defined as 

k nn (k) = Y / k'P(k'\k), (19) 

k> 

where the conditional probability P(k'\k) is the frequency that a node with 
degree k connects to a node with degree k! . Using Kronecker's delta function, 
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Fig. 5. Global clustering coefficient C with a = 5 and N = 12800. (A) 3D plot of 
the coefficient versus the variables, £ and m/a (axis direction changed for clarity). 
Projections of the plot to the (m/a)-C plane (B) and the £-C plane (C). 

we redefined the degree-degree correlation as 



~l2i=l{knn}i x $(ki k) 

NP(k) 



(20) 



where (k nn )i denotes the average nearest-neighbor degree, written as 

(k n n)i = T~ J]] (21) 
* heV(i) 



where V(i) corresponds to the set of neighbors of node i. 

Here, we give the numerical solutions for the degree-degree correlation of our 
model. Figure [6] shows the degree-degree correlations with a = 5, N = 50000, 
the different m, and £. The correlations follow the power law; this k nn (k) oc k u 
with — 1 < v < as a roughly observation. Negative and larger v tend to be 
seen for the larger m/a and £. 

Figure [6] (A) shows the degree-degree correlations with m = 2 and the different 
£. The exponent v decays with increasing £. Figure [6] (A) shows the degree- 
degree correlations with £ = 7 and the different m. The exponent v decays 
with increasing m. 
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Fig. 6. Degree correlations k nn (k) with a = 5 and ./V = 50000. (A) £ dependency 
with m = 2. A solid and dashed line correspond to oc fc~ 0,3 and oc k ' 2 , respectively. 
(B) m dependency with £ = 7. A solid line shows the relationship oc A;~ ' 8 . 



Due to the fitness updating and the PA, the degree-degree correlations follow 
the power law, reflecting the disassortativity of the networks. As previously re- 
ported, the disassortativity is not reproduced if we consider the PA mechanism 
only [27|I28] . 



8 Assortative coefficient 



The assortative coefficient (AC) [30] is can be thought of as a compendium 
parameter of the degree-degree correlations, and is defined as 

^ _ 4:{kikj) — (ki + kj} 2 



2(k 2 + k 2 ) - (ki + kj)' 



where ki and kj are the degrees of two vertices at the ends of an edge, and (• ■ ■) 
denotes the average over all edges. In other words, the AC is the correlation 
coefficient for the degree-degree correlation k nn (k), and takes values —1 < 
r < 1. The relationship between the AC and network structures is described 
as follows: 

i) In the case of r < 0, low-degree nodes tend to connect to high-degree 
nodes, indicating the disassortativity. Then, the degree-degree correlation 
k nn {k) decreases with increasing k. 

ii) In the case of r = 0, the degree-degree correlation k nn (k) is uncorrelated. 

iii) In the case of r > 0, high-degree nodes tend to connect to high-degree 
nodes, reflecting the assort at ivity. Then, the degree-degree correlation 
k nn (k) increases with degree k. 



12 



Fig. 7. Assortative coefficient r with a = 5 and N = 12800. (A) 3D plot of the 
coefficient versus the variables, £ and m/a. Projections of the plot to the (m/a)-r 
plane (B) and the £-r plane (C). Vertical axis of (A)-(C) is inverted for clarity. 

Figure [7J shows the numerical solutions of the AC for our model with a = 5, 
N = 12800, different m G [1,4], and £ e [0, 15]. 

Figure [7] (A) shows the AC with the different m/a and £ in three dimensions. 
The large negative AC generally tends to be seen for the larger m/a and £. 
Figure |7J (B) shows the projections of the plot to the (m/a)-r plane of Fig. 
[7J (A). Two biphasic curves (once r takes positive instead of negative) for 
the smaller £ 6 [0, 1] emerge because positive degree-degree correlations are 
exhibited in the case of m/a < 0.5 [27]. For larger £ G [3,15], assortative 
coefficients r monotonically decrease (upward) with increasing m/a. Figure [7J 
(C) shows the projections of the plot on the £-r plane for Fig. [7J (A). 



9 Comparison of statistical properties 

In order to validate our model we make a comparison of statistical proper- 
ties between biological networks and our model networks. We prepare two 
different types of networks; the gene regulatory [37] and the metabolic [S] 
networks of Escherichia coli. The gene regulatory network is represented as 
sets of graphs consisting of nodes and edges; they correspond to genes and 
interactions among genes, respectively. For simplicity, we extract the largest 
component from the networks, and replace directed and/or weighted edges 
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Fig. 8. Comparison of statistical properties between biological networks and model 
networks. The left column shows degree distributions P(k). Clustering spectra C(k) 
are in center column. Degree-degree correlations k nn (k) constitute the right column. 
Note that the degree-degree correlations k nn (k) are divided by (k 2 )/(k) for normal- 
ization. Symbols and dashed lines denote data of biological networks and our model 
networks, respectively. (I) Gene regulatory network in Escherichia coli [37] and our 
model with a = 3, m = 1, and £ = 20. (II) Metabolic networks in Escherichia coli 
[8] and our model with a = 3, m = 2, and £ = 1. The size of the model network is 
the same as total number of nodes of the biological network. 

with undirected and unweighted edges. Moreover, we remove multiple edges 
and self-loops. The metabolic network is transformed with the same proce- 
dures. 

In Fig. [8], we show the statistical properties of the biological networks and 
our model networks. We have good agreement with the real network's data 
and our model, demonstrating that our model reproduces the three striking 
statistical properties, also widely shared in the biological networks. In addition 
we confirm similar relationships with network data derived from the yeast [38] 
and the worm [39] (data not shown). 



10 Discussions 



Even though large amounts of structural data of biological networks have 
been available in recent years, little is known about the detailed mechanisms 
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behind the origin and evolution of the networks. However, we proposed a 
hypothetical model and its parameter space that can reproduce biological 
features from a statistical viewpoint. We think that the growth mechanism by 
merging modules and the FD-PA mechanism are important for evolutionary 
processes in biological networks. 

Since new genes commonly appear through complete or partial duplication 
of pre-existing genes [10114 1 f 42 f 43 1144] . the resultant sequences are subject to 
selective pressure as parts of the whole system and would give rise to changes in 
functional interactions among biomolecules, resulting in changes of statistical 
properties of network topology. 

In the case of metabolic networks, modules in our model can be regarded 
as a set of chemical compounds (nodes) and pathways (edges), and the merg- 
ing process may correspond to reminiscent tinkering processes such as bridging 
pathways between pre-existing nodes, enabling doubling pathways. These pro- 
cesses may positively or negatively contribute to increase robustness against 
mutations or give divergence and become a new pathway. The assumption of 
modular elements in our model is supported because module structures are 
found in many biological networks [23f4"5] . 

The original PA mechanism is exactly proportional to the number of edges 
of the node [16]. The FD-PA mechanism can consider additional or extra fit- 
ness of nodes. Rozenfeld et al. [M] and Zhu et al. [35] introduce into nodes 
the additional factors aging and a radius from a node in two-dimensional 
space. Because the additional factors are determined on the basis of a node's 
property, it is easy for us to infer that the effects from the factor is consid- 
ered to be underlie physicochemical properties of biomolecules such as mass, 
length, affinity, structural stabilities, and contents of bases or amino acids in a 
biomolecule. Together with the idea increasing fitness during the growth pro- 
cess reflects a certain variation of node's properties with time, and is relevant 
to specification and/or fixation process after duplication and divergence. 



11 Summary 

We have presented analytical and numerical results verifying that our model 
produces the three remarkable statistical properties that are widely shared in 
biological network structures. 

Scale-free connectivity — The degree distribution of our model follows power 
law, indicating the scale-free property; thus P(k) oc /c~ 7 . The degree exponent 
7 can be predicted by Eq. (fTTI) , and demonstrates the range 2 < 7 < 00 which 
includes the range of the degree exponents that are observed in biological 



15 



networks 2 < 7 < 3 [T0lll2|ll8j . The degree exponents take smaller values, 
meaning that the degree distribution becomes flat for the larger m/a and £. 

Hierarchical modularity — The clustering spectrum of our model follows the 
power law with the exponent — 1; hence C(k) oc k~ l , reflecting the hierarchical 
modularity [24.25J. The power-law spectrum is nearly independent of £ and 
m/a. The modularity of the whole network is higher for large £ and smaller 
m/a. The hierarchy manifests growing mechanism by merging modules. 

Disassortativity — The degree-degree correlation of our model follows the 
power law; therefore k nn (k) oc k v , reflecting the disassortativity. The exponent 
v has the range — 1 < v < which includes the range of the exponents that 
are seen in biological networks [T3"fl4[29| . The result indicates that the FD 
mechanism with positive £ is critical to reproduce this property. In addition, 
the disassortativity is sensitive to £ and m/a. 
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